## Authors: Alexander Herzog and Kenneth Benoit
## Date: May 30, 2015
## Replication file for JOP article "The Most Unkindest Cuts: Speaker Selection and Expressed Government Dissent During Economic Crisis"

rm(list = ls(all = TRUE))

library(rjags)
library(ggplot2)
library(texreg)

### PATHS ##############################
dataAll <- "./generated_data/working_data_all.RData"
dataSub <- "./generated_data/working_data_speakers.RData"
########################################

### DATA ##############################
load("./estimated_models/joint_model3_multilevel_posterior.RData")
########################################

# Note: Model 3 only

coef.names.sel <- c(
    "Intercept",
    "Constituency unemployment",
    "Electoral safety",
    "Legislative seniority",
    "Party leader",
    "Government",
    "Cabinet member",
    "Party size (log)",
    "Debate days (log)",
    "Crisis",
    "Const. unemployment x Crisis",
    "Electoral safety x Crisis",
    "Government x Crisis",
    "Const. unemployment x Government",
    "Electoral safety x Government",
    "Const. unemployment x Government x Crisis",
    "Electoral safety x Government x Crisis")

stats.names.sel <- c("$\\mu_{\\alpha}^{\\text{sel}}$", "$\\sigma_{\\alpha}^{\\text{sel}}$")

coef.names.out <- c(
    "Intercept",
    "Constituency unemployment",
    "Electoral safety",
    "Legislative seniority",
    "Government",
    "Cabinet member",
    "Crisis",
    "Const. unemployment x Crisis",
    "Electoral safety x Crisis",
    "Government x Crisis",
    "Const. unemployment x Government",
    "Electoral safety x Government",
    "Const. unemployment x Government x Crisis",
    "Electoral safety x Government x Crisis")

stats.names.out <- c("$\\mu_{\\alpha}^{\\text{out}}$", "$\\sigma_{\\alpha}^{\\text{out}}$", "$\\sigma$", "$\\rho$")


# Extract model coefficients and quantiles
# ----------------------------------------
model.results <- data.frame(
    variable = character(), 
    model = character(),
    coef = numeric(),
    ci.low = numeric(),
    ci.up = numeric())


s <- summary(joint3.posterior)

# Selection model
# ==============

# coefficients
# ------------
# intercepts and betas
as <- mean(s$statistics[grepl("as\\[",rownames(s$statistics)),][, "Mean"])
bs <- s$statistics[grepl("bs\\[",rownames(s$statistics)),][, "Mean"]

# mu.as
mu.as <- s$statistics[rownames(s$statistics)=="mu.as"][1]

# sigma.as
sigma.as <- s$statistics[rownames(s$statistics)=="sigma.as"][1]

# combine together
coef <- c(as, bs, mu.as, sigma.as)


# credible intervals
# ------------------
# intercepts and betas
ci.as <- colMeans(s$quantiles[grepl("as\\[",rownames(s$quantiles)),][, c("2.5%","97.5%")])

ci.bs <- s$quantiles[grepl("bs\\[",rownames(s$quantiles)),][, c("2.5%","97.5%")]

# mu.as
ci.mu.as <- s$quantiles[rownames(s$statistics)=="mu.as", c("2.5%","97.5%")]

# sigma.as
ci.sigma.as <- s$quantiles[rownames(s$statistics)=="sigma.as", c("2.5%","97.5%")]

# combine together    
ci <- rbind(ci.as, ci.bs, ci.mu.as, ci.sigma.as)  

# labels
variable <- c(coef.names.sel[1:(length(bs)+1)], stats.names.sel)

# add to data frame
d <- data.frame(model="selection", variable, coef, ci.low=ci[,1], ci.up=ci[,2])

model.results <- rbind(model.results, d)


# Outcome model
# =============

# coefficients
# ------------
# intercepts and betas
ao <- mean(s$statistics[grepl("ao\\[",rownames(s$statistics)),][, "Mean"])
bo <- s$statistics[grepl("bo\\[",rownames(s$statistics)),][, "Mean"]

# mu.as
mu.ao <- s$statistics[rownames(s$statistics)=="mu.ao"][1]

# sigma.ao
sigma.ao <- s$statistics[rownames(s$statistics)=="sigma.ao"][1]

# sigma
sigma <- s$statistics[rownames(s$statistics)=="sigma"][1]

# rho
rho <- s$statistics[rownames(s$statistics)=="rho"][1]

# combine together
coef <- c(ao, bo, mu.ao, sigma.ao, sigma, rho)


# credible intervals
# ------------------
# intercepts and betas
ci.ao <- colMeans(s$quantiles[grepl("ao\\[",rownames(s$quantiles)),][, c("2.5%","97.5%")])

ci.bo <- s$quantiles[grepl("bo\\[",rownames(s$quantiles)),][, c("2.5%","97.5%")]

# mu.ao
ci.mu.ao <- s$quantiles[rownames(s$statistics)=="mu.ao", c("2.5%","97.5%")]

# sigma.ao
ci.sigma.ao <- s$quantiles[rownames(s$statistics)=="sigma.ao", c("2.5%","97.5%")]

# sigma
ci.sigma <- s$quantiles[rownames(s$statistics)=="sigma", c("2.5%","97.5%")]

# rho
ci.rho <- s$quantiles[rownames(s$statistics)=="rho", c("2.5%","97.5%")]

# combine together    
ci <- rbind(ci.ao, ci.bo, ci.mu.ao, ci.sigma.ao, ci.sigma, ci.rho)  

# labels
variable <- c(coef.names.out[1:(length(bo)+1)], stats.names.out)

# add to data frame
d <- data.frame(model="outcome", variable, coef, ci.low=ci[,1], ci.up=ci[,2])

model.results <- rbind(model.results, d)


# distinguish between betas and stats
model.results$type <- ""
model.results$type[model.results$variable %in% c(stats.names.sel, stats.names.out)] <- "stat"
model.results$type[!(model.results$variable %in% c(stats.names.sel, stats.names.out))] <- "beta"


# factor ordering of variables for plots
n <- length(coef.names.sel)
model.results$variable <- factor(model.results$variable, levels=c(stats.names.out, stats.names.sel, coef.names.sel[order(seq(1,n,1), decreasing=TRUE)]))


# Plot coefficients
# -----------------
# selection model
d <- subset(model.results, model.results$model=="selection" & model.results$type=="beta")
d <- d[!(d$variable=="Intercept"),]
p <- ggplot(data=d) +
    theme_bw() +
    geom_pointrange(
        aes(x=variable, y=coef, ymin=ci.low, ymax=ci.up),
        lwd = 1/2, position = position_dodge(width = 1/2),
        shape = 21, fill = "black") +
    geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
    coord_flip() +
    xlab("") +
    ylab("Coefficients and 95% credible intervals") +
    theme(text = element_text(size=16)) +
    theme(axis.title.x = element_text(size=12, vjust=-0.5)) +
    ggtitle("Speaker selection")

pdf("./plots/coefplot_joint3_selection.pdf", 8, 6)
print(p)
dev.off()


# outcome model
d <- subset(model.results, model.results$model=="outcome" & model.results$type=="beta")
d <- d[!(d$variable=="Intercept"),]
p <- ggplot(data=d) +
    theme_bw() +
    geom_pointrange(
        aes(x=variable, y=coef, ymin=ci.low, ymax=ci.up),
        lwd = 1/2, position = position_dodge(width = 1/2),
        shape = 21, fill = "black") +
    geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
    coord_flip() +
    xlab("") +
    ylab("Coefficients and 95% credible intervals") +
    theme(text = element_text(size=16)) +
    theme(axis.title.x = element_text(size=12, vjust=-0.5)) +
    ggtitle("Position taking")

pdf("./plots/coefplot_joint3_outcome.pdf", 8, 6)
print(p)
dev.off()


# Generate table
# --------------
load(dataAll)
load(dataSub)

N <- nrow(dataAll)
N.groups <- length(unique(dataAll$m))
N.speakers <- length(unique(dataSub$m))


sel <- with(model.results[model.results$model=="selection",],
            createTexreg(
                coef.names=as.character(variable),
                coef=coef,
                ci.low=ci.low,
                ci.up=ci.up,
                gof=c(N, N.groups, N.speakers),
                gof.names=c("N obs.", "N legislators", "N speakers"),
                gof.decimal=c(FALSE, FALSE, FALSE))
            )

out <- with(model.results[model.results$model=="outcome",],
            createTexreg(
                coef.names=as.character(variable),
                coef=coef,
                ci.low=ci.low,
                ci.up=ci.up,
                gof=c(N, N.groups, N.speakers),
                gof.names=c("N obs.", "N legislators", "N speakers"),
                gof.decimal=c(FALSE, FALSE, FALSE))
            )


texreg(
    list(sel, out),
    custom.model.names = c("Speaker selection","Position taking"),
    #reorder.coef = c(1,2,3,4,5,6,7,8,11,12,13,14,15,16,17,9,10),
    digits=2,
    single.row=FALSE,
    leading.zero=TRUE,
    caption.above=TRUE,
    booktabs=TRUE,
    bold=0.05,
    use.packages=FALSE,
    table=FALSE,
    stars=0,
    ci.test = NULL,
    ci.force=TRUE,
    center=TRUE,
    file="./tables/table_joint3_model.tex"
)

